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1.  Introduction 

Suppose  we  have  two  stochastic  systems  (perhaps  alternative  designs 
for  a new  system)  which  are  to  be  compared.  Assume  that  these  systems 
are  represented  by  two  regenerative  processes  X(  i)  = (X  ( i)  : t > 0) 
for  i = 1,  2;  see  CRANE  and  LEMOINE  ( I977)  or  IGLEHART  ( I978)  for  a 
discussion  of  regenerative  processes  and  their  role  in  simulation.  Under 
mild  regularity  conditions  the  distribution  of  X^( i)  converges  to  the 
distribution  of  some  limiting  random  variable  (or  vector)  X( i) ; this  type 
of  convergence  is  known  as  weak  convergence  and  written  X^(  i)  =>  X(  i) 
as  t t oe.  Simulators  often  speak  of  X(  i)  as  the  steady-state  configura- 
tion of  system  i and  take  as  the  performance  criterion  of  the  system 
tf  = E(f^(X(i))},  where  f^.  (i  = 1,2)  are  given  real-valued  functions 
defined  on  the  state  space,  E(i),  of  process  X(  i) . When  comparing  the 
two  systems,  we  wish  to  estimate  the  sign  of  r^  - by  constructing  a 
confidence  Interval  for  the  quantity. 

This  research  was  supported  by  National  Science  Foundation  Grant 
MCS-23607  and  Office  of  Naval  Research  Contrac;t  NOOOUt-76-C-(^78. 


Of? 


Simulation  folklore  suggests  that  "common  random  numbers"  be  used 


in  this  situation  in  order  to  reduce  variance;  see  FISHMAN  (1975),  Section 
11.7^  and  KLEIJNEN  ( l97k)  for  a discussion  of  this  technique.  The  basic 
idea  here  is  to  induce  a positive  correlation  between  the  two  systems  by 
simulating  the  two  systems  with  a common  sequence  of  random  numbers.  Not 
only  does  one  save  the  computer  time  required  to  generate  the  second 
sequence  of  random  numbers  (which  would  be  required  if  the  two  systems 
were  simulated  independently)^  but  the  confidence  interval  for  r^^  - r^ 
will  be  shorter  provided  the  positive  correlation  mentioned  above  is 
achieved.  Inspite  of  the  generally  knowledged  appeal  of  this  technique, 
we  know  of  few  published  studies  which  actually  carry  out  the  technique 
and  document  the  savings  to  be  expected;  one  such  study  is  that  of 
MITCHELL  (1975). 

When  the  processes  being  compared,  X( 1)  and  X(2),  are  regenerative, 

we  are  able  to  provide  a rigorous  (asymptotic)  analysis  of  the  comparison 

technique  described  above.  This  we  do  in  Section  2 of  the  paper.  While 

positive  correlation  is  normally  expected  when  using  common  random  numbers, 

nothing  is  guaranteed  in  this  respect.  One  always  enjoys  the  economy  of 

having  to  generate  only  half  as  many  random  numbers,  however,  the  variance 

reduction  achieved  may  be  minimal  or  even  a variance  addition.  Conditions 

for  obtaining  the  desired  positive  correlation  are  discussed  in  Section 

Section  k is  devoted  to  several  simple  stochastic  systems  which  were 

actually  simulated.  Here  we  are  able  to  see  exactly  what  common  random 

numbers  buys  the  simulator  in  terms  of  increased  computational  efficiency. 

Finally,  in  Section  ^ we  state  our  basic  conclusions  about  the  use  of  common 
% 

random  numbers  in  comparing  stochastic  systems  using  regenerative  simulation. 
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2.  Comparing  Regenerative  Processes 


Let  X( 1)  and  X(2)  be  the  two  regenerative  processes  introduced 
in  Section  1.  Assume  X^(  i)  ^ X(  i)  as  t T <»  for  i = 1^2.  Given 
: E(  i)  -*R,  we  wish  to  construct  a confidence  interval  for  r^^  - r^. 

The  regenerative  processes  arising  most  commonly  in  simulations  are 
discrete  and  continuous  time  Markov  chains  and  semi-Markov  processes  all 
of  which  are  positive  recurrent.  An  efficient  method  for  reducing  the 
simulation  of  a continuous  time  Markov  chain  (M.C.)  to  the  simulation  of 
a related  discrete  time  M.C.  was  presented  in  HORDIJK,  IGI£HART,  and 
SCHASSBERGER  ( I976) . The  same  method  can  be  applied  to  semi-Markov 
processes;  see  IGLEHART  ( I978) . In  our  simulations  presented  in  Section  k 
this  method  is  used.  Hence  to  focus  our  attention  on  the  comparison 
problem  at  hand^  we  assume  that  both  X( 1)  and  X(2)  are  irreducible^ 
aperiodic^  positive  recurrent  Markov  chains  in  discrete  time.  Under 
these  conditions  X^( i)  =»  X( i)  as  n t where  X(  i)  has  the  stationary 
(and  steady-state)  distribution  i)  = (rr^Ci)  : j € That  is, 

P(X(i)  = j)  = itj(i)  for  jeE(i).  Then 


r = E{f  (X(i))}  = E fjj)  Tt  (i)  , i = 1,2  . 

jeE^  ^ J 
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Assume  for  this  discussion  that  E = E( 1)  = E(2)  and  that  the  state 
0 C E;  this  is  no  restriction  on  our  method  only  a notational  convenience. 
Then  let  Xq( i)  = 0,  Tq( i)  = 0,  and  define  the  m^^  entrance  to  state  0 
by  X(i)  to  be 


The  eymbol  E^t*)  is  short  for  the  conditional  expectation  E(-|XQ(i)  = 0}. 
For  more  background  on  the  regenerative  method  see  CRANE  and  LEMOINE 
( 1977)  IGLEHART  ( I978) . Let  = EQ{Z^(i))  which  we  assume  is 
positive  and  finite.  Then  two  central  limit  theorems  (c.l.t.'s)  follow 
from  the  regenerative  structure  of  these  positive  recurrent  Markov  chains. 


k 


One  Is  based  on  Che  number  of  0-cycles  simulated  and  Che  other  on  Che 
number  of  steps  of  the  M.C.  simulated.  We  define  two  point  estimates 


of  r^  by 


Wi)l  , 


where  n is  the  number  of  0-cycles  simulated  and  N is  the  number  of 
steps  of  the  M.c.  simulated.  The  two  c.l.t.’s  are  the  following: 


(2.1)  n^  ^[r^(n)  - r^]/(  a^/EQ(T^(  i) ))  N(0, 1) 


(2.2)  N^'^2rr^(N)  - r^]/(a^/Ey^{T^(i)))  =>N(0,1) 


as  n and  N t «.  Either  (2.1)  or  (2.2)  can  be  used  to  construct  a 
confidence  interval  for  r^. 

Now  suppose  we  wish  to  construct  a confidence  interval  for 
by  simulating  the  two  processes  X( 1)  and  X(2)  independently;  i.e., 
independent  sequences  of  random  numbers  will  be  used  to  generate  the 
sample  paths  of  the  two  processes.  Form  the  vectors  r = (r^^,  rg)  and 
r(N)  = (rj^(N),  rg(N)).  Then  we  can  easily  obtain  from  (2.2)  the  bivariate 


c.  1. c. 


(2.3) 


N^''^[£(N)  - r]  =>N(0,  A)  , 
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where  N(0,  A)  is  a two-dimensional  normal  vector  with  mean  vector 
0 = (0,0)  and  covariance  matrix 

/ 0 

0 

An  application  of  the  continuous  mapping  theorem  [BILLINGSLEY  (I968), 
Theorem  5.1]  to  (2.5)  yields  the  following  c.l.t.  which  can  be  used  to 
construct  a confidence  interval  for 

(2.4)  N^/2[(r^(N)  - r2(N))  - (r^  - r^)  ]/o  =»  N(0, 1)  , 

where 

2 2 
2 ‘'l  ‘^2 

^ = Eq{t^(1))  ^ Eq{t^(2))  • 

A c.l.t.  comparable  to  (2.4)  but  based  on  0-cycles  can  also  be  obtained. 

Our  goal  in  using  common  random  numbers  to  generate  the  sample 
paths  of  X(  1)  and  X(2)  is  to  produce  a shorter  confidence  interval 
for  length  of  simulation  run  (number  of  steps  of 

the  Markov  chains  generated).  In  other  words  we  seek  a c.l.t.  similar 
to  (2.4)  but  with  a smaller  value  of  a.  To  accomplish  this  we  generate 
the  bivariate  M.c.  X = (X  : n > O),  where  X = (X  ( 1) , X (2) ) . At 
each  jump  of  the  process  X the  same  random  number  is  used  to  generate 
the  Jumps  of  the  two  marginal  chains  X(  1)  and  X(2).  The  marginals 
of  the  process  X are  seen  to  have  the  same  finite-dimensional 
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distributions  as  the  original  chains  X( 1)  and  X(2);  however^  the 
marginal  chains  are  now  dependent.  The  state  space  of  the  chain  X 
is  denoted  by  F which  is  a (possibly  proper)  subset  of  E x E.  We 
assume  here  that  the  chain  X is  also  irreducible,  aperiodic,  and  positive 
recurrent.  These  conditions  are  not  automatic  but  will  usually  hold  for 
practical  simulations.  Furthermore,  we  assume  for  convenience  that 
(0,0)  € F and  use  Chat  state  to  form  regenerative  cycles.  Mote  that 
X^  X as  n 00  and  Che  marginal  distributions  of  X are  the  same  as 

those  of  X(  1)  and  X(2),  namely,  {jtj(i)  : j € E)  for  i = 1,2.  For 
any  real-valued  function  f : F -♦  R satisfying  E(|f(X)(]  < o»  the 
regenerative  method  can  be  applied  to  X to  estimate  E(f(X)}.  Let 
Xq  = (0,0)  and  form  (0,0)-cycles  which  begin  at  the  times  T^  = 0, 


T = inf(n  > T , : X = (O.O))  . m > 1. 

m m- 1 n ' ^ ^ ' — 


Also 


th 


let  T = T - T ,.m>l,  be  the  length  of  the  m cycle  and 
m m m- 1'  — ' 


n=l  , 
m- 1 


m > 1 . 


Here  the  f-functions  are  f(j,k)  = fj^(j)  and  f(j,k)  = f2(k).  Set 

Z'(i)  3 Y'(i)  - r,  T . Since  the  ratio  formula  still  holds  for  Che 
m'  ' m ' i m 

process  X,  (Z^( i) ) =0  for  i=  1,2.  Let 


"ij  " . 


I,j  . 1,2 


i ; 
1 


which  we  assume  is  finite  and  non-zero.  Since  the  vectors 

Z'  = (Z'(l).  Z'(2))  are  i.i.d.,  the  standard  c.l.t,  yields 
~m  m ' ^ m'  ' ' ’ 


(2.5) 


n‘  Z Z ' N(  0,  E)  , 


m=l 


where  E = Just  as  we  are  able  to  go  from  (2,1)  to  (2,2)  in 

the  one -dimensional  case  it  is  possible  to  obtain  from  (2.5)  the  c.l.t. 


(2.6) 


N^'^^[r(N)  - r]  N( 0, B)  , 


where  B = argument  leading  to  (2,6)  is  essentially 

the  same  as  that  given  by  CHUNG  (i960),  Theorem  I6.I.  Again  using  the 
continuous  mapping  theorem  in  conjunction  with  (2.6)  yields 


(2.7) 

where 


N^'^i;(r^(N)  - rgCN))  - ( r^-r^)  ]/v  =>  N(0, 1)  , 


V - (<^11  + <^22  ■ ^‘^12^/^{0,0)^’^1^  • 


A c.l.t.  comparable  to  (2.7)  but  in  terms  of  n (0,0) -cycles  of  X can 
also  be  obtained.  Now  consider  the  marginals  of  (2.6)  in  conjunction 
with  (2.2).  Since  the  marginals  of  the  chain  X have  the  same  stochastic 
structure  as  the  chains  X(  1)  and  X(2)  considered  separately,  these 
two  c.l.t. 's  must  be  identical.  Hence 
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2 2 

Thus  upon  comparing  the  constant  a in  (2.1).)  and  v in  (2.7)  we 
2 2 

conclude  that  v < a if  and  only  if  Oj^g  ^ Section  3 we  will 

examine  conditions  on  the  functions  f and  processes  X( i)  which 
guarantee  that  <^12  ^ 

The  measure  of  variance  reduction  we  use  is 


SOj  for  examole,  if  R = 0.5^  then  only  half  as  many  steps  of  the  Markov 
chain  X need  be  simulated  to  obtain  a confidence  interval  of  specified 
length  for  would  be  required  when  simulating  X(  1)  and  X(2) 

independently.  In  addition^  of  course,  only  one  stream  of  random  numbers 
need  be  generated.  While  we  have  worked  here  with  discrete  time  Markov 
chains,  the  same  method  can  be  used  for  continuous  time  Markov  chains, 
semi-Markov  processes,  and  discrete  time  Markov  processes  with  a general 
state  space.  The  examples  treated  in  Section  k illustrate  the  effectiveness 
of  the  method  when  applied  to  a variety  of  these  stochastic  processes. 


r 

I 
r 

! 3.  Guaranteeing  Variance  Reductions 

In  this  section  we  investigate  conditions  under  which  the  variance 
reduction  obtained  when  > 0 can  be  guaranteed.  Our  major  result  is 

that  if  and  fg  are  monotonic  functions  (in  the  same  direction) 

and  if  2(1)  and  ^(2)  satisfy  a stochastic  monotonicity  condition^ 
then  Oj^  > 0.  This  result  is  related  to  other  work  on  monotonicity 
and  antithetic  variates;  (see  ANDREASSON  (1972),  MITCHELL  (1973),  and 
KLEIJNEN  ( 19714-)).  When  using  antithetics  (if  U is  a random  variable 
uniformly  distributed  on  [0,1]  then  U and  1-U  are  said  to  be  an 
antithetic  pair)  one  is  generally  interested  in  only  one  stochastic 
process,  not  in  comparing  the  output  of  two  or  more  processes.  Also  in 
the  antithetic  scheme  variance  reductions  are  obtained  by  creating  negative 
correlation  rather  than  the  positive  correlation  we  seek  here.  If  only 
one  process  is  to  be  considered,  the  sample  paths  of  ^(2) 

may  be  generated  using  antithetic  streams  of  random  numbers.  Under 
proper  conditions  the  results  of  this  section  may  then  be  applied  to 
guarantee  the  desired  negative  correlation  ( provided  that  the  two  dimen- 
sional process  X Is  regenerative). 

The  notion  of  associated  random  variables  can  be  used  to  guarantee 
nonnegative  correlation.  The  following  definition  and  properties  may  be 
found  in  ESARY,  PROSCHAN  and  WALKUP  ( I967) . 

(3.1)  DEFINITION.  Random  variables  T = (Tj^,  ...,  T^)  are  said  to  be 
associated  if 

cov{f(T),  g(T)  } > 0 , 
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for  all  nondecreasin^  functions  f and  g for  which  E[f(T)}^  E[g(T)), 
and  E{f(T)  g(T))  exist. 


(3.2)  PROPERTY.  Any  subset  of  associated  random  variables  are  associated. 

(3.3)  PROPERTY.  If  two  sets  of  associated  random  variables  are  independent 
of  one  another,  then  their  union  is  a set  of  associated  random  variables. 

(3.ii-)  PROPERTY.  The  set  consisting  of  a single  random  variable  is 
associated. 

(3.5)  PROPERTY.  Nondecreasing  functions  of  associated  random  variables 
are  associated. 

A class  of  processes  for  which  nonnegative  correlation  can  be 
guaranteed  is  stochastically  monotone  Markov  chains  (s.m.m.c.).  This 
class  was  introduced  by  DALEY  ( I968)  and  includes  many  of  the  basic 
queueing  models  such  as  the  waiting  time  process  in  the  GI/G/l  queue 
and  the  embedded  Markov  chains  used  to  study  the  M/G/1  and  g/m/s 
queues.  In  the  following  definition  let  i be  a fixed  index. 

(3.6)  DEFINITION.  Let  X( i)  = (X^(i),  n > 0}  be  a real  valued  Markov 

process  with  initial  distribution  Pj^(x)  = P{XQ(i)  < x)  and  transition 
function  Pj^(x,A)  = *■)  £ A!x^(i)  = x)  (for  measureable  sets  A). 

X(  i)  is  said  to  be  a stochastically  Monotone  Markov  chain  if  for  every 
y,  (-“,y])  is  a nonincreasing  function  of  x. 
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Define  the  inverse  distribution  functions 


P"  (u)  = inf{y  : P^(y)  > u) 


(x,“)  = inf(y  : P^(x,  ( -»,  y])  > u}  . 


Notice  that  if  X(  i)  is  a s.tn.m.c.  then  P^  (x,u)  is  an  increasing 
function  in  both  arguments.  This  fact  will  enable  us  to  show  that  for 
each  n>0  {Xq( 1) , . . . , Xj 1) , Xq( 2) , . . . , XJ 2) } are  associated. 

We  shall  henceforth  assume  that  the  sample  paths  of  X(  i)  are 
generated  on  the  computer  using  the  inverse  transformation  scheme 


C3.7)  XqU)  = p-\Uq)  , 


(5.8) 


X (i)  = P'^X  ,(i),  U ) , n > 1 

n'  ^ i ' n- 1'  ' » n'  ’ — 


where  U = {U^,  n > 0)  is  a sequence  of  pseudorandom  numbers.  U is^ 
of  course^  assumed  to  be  a sequence  of  i.i.d.  random  variables  uniformly 
distributed  on  [0,1], 


(3.9)  THEOREM.  If  X( 1)  and  X( 2)  are  both  stochastically  monotone 
Markov  chains  with  sample  paths  generated  by  (3.7)  and  (3.8),  then  for 
each  n >0  (Xq(1),  ...,  X^( 1) , Xq(2),  ...,  X^(2))  are  associated  random 

variables. 


P 


PROOF.  The  proof  is  by  inducclon.  For  n = 0 (^.h)  implies  chat  [U^) 

is  associated  and  since  ^ nondecreasing  function  of  for 

each  i,  (5.5)  yields  that  {Xq{1),  ^^(2))  are  associated.  Assume  now 
that  {Xq(1)^  ...^  associated.  Since 

is  independent  of  this  set,  (Xq(1),...  , \(  1) , ' • • • ' » ^n+l^ 

are  associated  (by  (5.5)).  The  map  which  Cakes  these  random  variables 
into  {Xq(1),  ...,  X^(l),  X^^^(l),  Xq(2),  ...,  X^(2),  X^^^(2)}  is  non- 
decreasing  because  X( 1)  and  X(2)  are  both  s.m.m.c. 's.  Property  (5-5) 
Chen  yields  Che  final  result.  □ 

The  previous  theorem  can  now  be  used  to  show  that  when  simulating 
s.m.m.c. 's  using  common  random  numbers  a reduction  in  variance  is  obtained. 


(5.10)  THEOREM.  Let  X( 1)  and  X( 2)  both  be  stochastically  monotone 
Markov  chains  with  sample  paths  generated  by  (5*7)  ®tid  (5.8).  Let  fj^ 
and  be  nondecreasing  functions.  If 


(i) 


< 00 


(ii) 


e{( 


< 00 


for  i = 1,  2, 


then  a^2  ^ 0* 
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PROOF.  Let 


n 

S ( i)  = Z S ( i)  is  then  a nondecreasing 

n Ik  n 

function  of  associated  random  variables  so  that 

associated.  Therefore  covfS  (1),  > 0.  (This  covariance  exists 

and  is  finite  by  ( i)  and  (ii),  see  SMITH  (1955)-)  Theorem  8 of  SMITH 
(I955),  which  may  be  applied  under  assumptions  ( i)  and  (ii),  implies  that 

lim  i covfSjl),  SJ2)}  = 0^2  , 
n — » 00 


so  that  aj^2  ^0* 

If  f and  fg  are  both  decreasing  then  Oj^g  ^ 
and  fg  are  monotonic  in  opposite  directions,  that  is  if  one  is  increasing 
and  the  other  is  decreasing,  then  ^ case  antithetics 

should  be  used  to  ensure  a^g  > 0) . 
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U.  Examples 


In  this  section  we  investigate  the  magnitude  of  variance  reductions 
obtained  when  using  common  random  numbers  in  simulations  of  some  simple 
stochastic  models.  These  models  are  the  waiting  time  process  in  the 
M/G/I  queue  and  the  queue  length  processes  in  the  finite  capacity  M/m/s 
queue,  repairman  problem,  and  a repairman  problem  with  two  repair  facilities 
(sometimes  called  the  central  server  model).  For  these  models  the  method 
was  able  to  produce  quite  substantial  variance  reductions.  In  addition 
we  investigated  the  use  of  common  random  numbers  when  comparing  two 
different  (3,S)  inventory  policies.  To  our  surprise  variance  reductions 
were  slight  in  this  case. 

The  Single  Server  Queue 

Let  i)  waiting  time  of  the  nth  customer  in  the  ith 

GI/G/l  queue  which  we  wish  to  study.  Let  {S^(i)  : n > 0)  be  the 
sequence  of  i.i.d.  service  times  (with  mean  and  distribution 

function  G^)  and  {A^( i)  : n > 1)  be  the  i.i.d.  interarrival  times 
(with  mean  and  distribution  function  F^)  for  this  queue.  Set 

X^(  i)  = j^(  i)  - A^(  i) . The  waiting  times  are  then  defined  by 

Wo( i)  = 0 

W ,(i)  = [W  (i)  + X .(i)!"^  , n > 0 

n+1'  ''  n'  ''  n+1^  ' ' - 

where  for  any  real  number  a,  [a)'*'  denotes  the  maximum  of  0 and  a. 

* Let  If  < 1 then  W^(  i)  W(  i)  • We  shall  be  interested 

I 

i 
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r 


in  escimating  E(W(1)}  - E(W(2)).  Recall  that  E{W( i) ) is  finite  if  j 

E{X^(i)  } < 00,  see  KIEFER  and  WOLFOWITZ  ( I956)  . Let  {V^  : n > 0)  and  | 

(U^  : n > 1)  be  independent  sequences  of  i.i.d.  uniformly  distributed  I 

random  variables  which  generate  the  service  and  interarrival  times  by 


s.d) 

= «'l  (’n)  - 

n > 0 

A„(l) 

■ . 

where 

C''  and  F"^ 

are  the  inverse 

distribution  functions  of 

and  F 

^ respectively. 

The  following 

theorem  states  conditions  under 

which  the  two  dimensional  process  W = {(W^(l),  W^(2)),  n > 0)  will  be 
regenerative. 

(L.JI  THEOREM.  Let  o.  < 1 and  E{X^(i)^)  <00  for  i = 1,2.  If  the 

1 n j 

joint  distribution  of  (^^^(l),  \(^))  ^as  a positive  density  in  an  open  | 

neighborhood  of  (0,0),  then  W Is  regenerative. 

i 

PROOF.  Let  e > 0 be  given  and  consider  the  discretized  waiting  time  ! 

processes 

i 

) 

“o< ° - I 

where  X^( i)  = ke  if  (k-l)e  < X^( 1)  < for  some  Integer  k.  The  i 

process  = {(W^(l),  W^(2)),  n > 0}  is  a Markov  chain  with  a countable 
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so  that  returns  to 


state  space.  Since  X ( i)  < X'^(  i) , W ( i)  <W^(i) 

(0,0)  occur  more  frequently  for  W than  for  W~,  The  condition  on 

the  joint  distribution  of  ^j,(  1)  ^n^^^  ensures  that  for  small 

enough  €,  will  be  irreducible  and  aperiodic  (so  that 

jr^(i,j)  = lim  P(W^(  1)  = ie,  W^(2)  = je)  exists).  It  therefore 
n — ► 00 

sufficies  to  show  chat  is  positive  recurrent. 

Let  T^  be  the  mth  time  W"  enters  (0.0).  We  seek  to  show 
that  E(T^)  < 00  and  since  jt^(0,0)  = l/E(Tp  we  need  only  show 
jT^(0,0)  > 0.  Let  e be  chosen  small  enough  so  that  the  traffic 
intensity,  p^,  in  each  discretized  queue  is  less  than  one  and  so  Chat 
EC(X^(  i)"^)^}  < 00.  Then  W^(i)=>W'(i)  and  E(W^(i))<oo.  Since  < 1, 

(l^.L)  0 <Tr^(0)  = Urn  P{W^(  1)  = 0} 

n 00 

00 

= lim  Z P{W^(1)  = 0,  W^(2)  = kc)  . 
n -*  00  k=0 

Then 

P(W^  = (0,  ke)}  < P{W^(2)  = ke} 

< P(W^(2)  > k€) 

<P{W"(2)>k€)  , 

the  last  inequality  being  true  because  the  Gl/G/1  queue  is  a s.m.m.c. 

00 

But  y,  P(W^(2)  > k€ ) <00  (since  E{W'(2))  < oo)  so  that  by  the  dominated 
k=0 

convergence  theorem  we  may  interchange  limits  and  summation  in  (k.k). 
Therefore 
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I 


I 


0 < 7,  n^(0,k) 
k=0 

and  since  n^(i,j)  is  either  0 for  all  i and  j or  greater  than 
0 for  all  i and  j it  must  be  true  that  jr^(O^O)  >0.  □ 

This  theorem  may  be  extended  to  multiple  server  queues  and  to 
situations  in  which  more  than  two  queues  are  being  considered.  In 
addition  these  queues  possess  the  proper  monotonicity  characteristics 
to  ensure  that  > 0. 

The  first  two  sections  of  Table  1 report  variance  reductions  when 
two  M,''G/1  queues  are  compared.  The  service  times  were  chosen  to  have  the 
Weibull  distribution: 

G^(x)  = 1 - exp(-(r^x)°'^)  , x>0  I 

for  constants  > 0.  If  = 1,  G^  reduces  to  the  exponential 

distribution.  The  figures  in  Table  1 are  point  estimates  and  90^  con- 
fidence intervals  based  on  N Independent  replications  of  C cycles  of 

the  two  dimensional  process  W.  For  example,  in  the  section  of  Table  1 

2 2 
we  estimate  R to  be  .11^  (and  a 90^^  confidence  interval  for  R is  j 

1 

.II4.2  J-  .015).  Notice  that  cycles  for  W are  not  much  longer  than  those  | 

of  the  individual  processes.  In  fact  for  the  processes  compared  in  the  1 

first  section  of  Table  1 it  can  be  shown  that  W (1)  < W (2)  for  all  n 

n'  ' — n 

so  that  T = T (2).  The  random  number  generator  described  in  LEARMONTH 
m m'  ' 

and  LEWIS  (197^)  was  used  for  all  simulations  reported  in  this  paper. 
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Point  Estimates  and  90‘j^  Confidence  Intervals  Using 


Continuous  Time  Markov  Chains 


The  use  of  common  random  numbers  in  comparing  two  or  more  continuous 
time  processes  is  limited  by  problems  in  the  "synchronization"  of  the 
random  number  streams  (see  KLEIJNEN  (197^)*  This  problem  can  be  overcome 
in  the  case  of  continuous  time  Markov  chains  and  semi-Markov  processes 
by  transforming  the  continuous  time  processes  into  appropriate  discrete 
time  Markov  chains.  Details  of  this  transformation  are  given  in  HORDIJK, 
IGLEHART  and  SCHASSBERGER  ( 197^) • Once  in  discrete  time  common  random 
numbers  may  be  used  to  generate  the  sample  paths  of  the  two  processes. 

This  procedure  has  been  used  to  investigate  variance  reductions  for 
three  finite  state  space  continuous  time  Markov  chains.  Because  the 
state  spaces  are  finite  the  multidimensional  processes  will  always  be 
positive  recurrent  (assuming  irreducibility) . 

The  first  two  examples^  the  queue  length  processes  in  the  finite 
capacity  M/m/s  queue  and  the  repairman  problem  with  spares  are  both 

city  M/m/s  queue  has  birth 

0 < i < M 

i > M 

1 < i < s 
s < i < M , 


birth  and  death  processes.  The  finite  capa 
and  death  parameters 


\ = 


^^i  = 


sn  , 


20 


where  M is  the  capacity  of  the  queue.  For  this  model  let  p = 'K/s\i, 
The  repairman  problem  has  parameters 


0 < i < m 

m < i < nn-n 

1 < i < s 

s < i < tin-n 

where  n is  the  number  of  operating  units^  m is  the  number  of  spare 
units^  s is  the  number  of  repairmen  and  \ and  ^ are  the  failure 
and  repair  rates  respectively  of  the  units.  Calculated  variance  reduc- 
tions for  these  two  models  are  reported  in  Tables  2 and  Ic  should 
be  noted  that  for  our  choice  of  parameters  1)  ^ n > 0 

provided  that  Xq(  1)  = Xq(2)  = 0.  These  are  examples  in  which  the  two 
dimensional  process  is  not  irreducible.  In  such  cases  attention  must  be 
focused  on  only  one  irreducible  class  of  states. 

The  next  example  is  a multidimensional  repairman  problem.  This 
example  can  be  modelled  by  the  closed  queueing  network  pictured  in 
Figure  1.  A more  detailed  description  of  this  model  may  be  found  in 
IGLEHART  and  LEMOINE  ( 197*^)  • The  parameters  chosen  for  this  model  were 
n = 10,  m a k,  X = 1,  p = .2,  s^  = 2,  = 5,  = 1.  The  effect  of 

varying  s^  on  the  mean  number  of  failed  units  was  studied.  Since  each 
individual  process  has  a two  dimensional  state  space,  a four  dimensional 


\ = 


nX  , 

(n+m-i)X  , 

su  , 


1 
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TABLE  2 

Calculated  Variance  Reductions  for  Comparing  Two  Finite  Capacity 
M/M/s  Queues^  Capacity  M = 15 
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TABLE  3 

Calculated  Variance  Reductions  for  Comparing 
Two  Repairman  Problems,  n = 10,  m = U,  X = 1. 
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FIGURE  1 


Repairman  Problem  with  Two  Repair 
Facilities  and  Spares 


2k 


state  space  is  obtained  when  using  common  random  numbers.  The  third 
section  of  Table  1 reports  the  estimated  variance  reduction  for  this 


example.  The  expected  cycle  length  for  the  four  dimensional  process  is 
much  larger  than  for  either  of  the  two  dimensional  processes^  but  it  is 
not  as  long  as  we  had  originally  feared.  Of  course  the  expected  cycle 
length  is  a function  of  the  return  states  chosen  and  care  must  be  taken 
so  that  regenerations  do  not  occur  too  infrequently. 


Inventory  Policies 

Our  final  example  is  comparing  two  different  (s^S)  inventory 
policies.  The  use  of  common  random  numbers  should  be  particularly 
well  suited  to  inventory  problems  since  it  intuitively  seems  better 
(i.e.j  less  variable)  to  compare  different  policies  by  subjecting  them 
to  the  same^  rather  than  independent^  demand  processes.  However^  the 
figures  in  Table  U indicated  that  very  little  variance  reduction  is 
obtained  for  this  model.  This  apparently  is  due  to  the  fact  that  (s^S) 
policies  grossly  violate  the  monotonicity  conditions  of  the  previous 
section. 

Let  denote  the  level  of  stock  at  the  beginning  of  the  nth 

period  and  let  be  the  demand  during  the  nth  period.  Then 


n+ 


l(i)  = 


X (1)  - D . 
n'  ' n ' 


if  X (i)  - D < s , 
n'  ' n > 

if  X (i)  - D < s . 
n'  ' n — 
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TABLE  4 

Calculated  Variance  Reductions  for  Comparing  Two  j 

(s^S)  Inventory  Policies,  r.  = E[X(i)]  | 

I 
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We  choose  demands  to  have  either  a geometric  or  Poisson  distribution  V7ith 
parameters  .5  and  2 respectively. 

The  discontinuity  at  s prevents  us  from  obtaining  strong  positive 
correlation  (and  hence  good  variance  reductions).  If  1)  near 
Sj  but  is  well  above  Sg  then,  with  high  probability,  the  first 

process  will  increase  while  the  second  will  decrease.  This  tends  to 
create  negative  correlation  rather  than  the  desired  positive  correlation. 

We  were  able  to  increase  the  variance  reductions  somewhat  by  using  antithetics 
when  in  certain  regions  of  the  two  dimensional  state  space.  However,  the 
generality  of  such  a procedure  seems  limited  since  it  may  be  quite  difficult 
and  costly  to  identify  regions  of  this  type  which  would  lead  to  good 
variance  reductions. 
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5 . Cone luslons 

In  this  paper  we  have  shown  how  the  method  of  common  random  numbers 
may  be  used  in  certain  regenerative  simulations  to  obtain  variance  reduc- 
tions. In  some  cases  substantial  variance  reductions  have  been  obtained, 
but  it  seems  reasonable  to  expect  that  as  the  complexity  of  the  processes 
being  simulated  increases  the  amount  of  variance  reduction  will  decrease. 
It  is  anticipated  that  the  primary  difficulty  in  the  implementation  of 
this  method  will  be  the  relatively  long  expected  cycle  length  for  the  two 
dimensional  process  X.  Since  the  validity  of  the  confidence  intervals 
formed  will  in  general  depend  upon  the  number  of  cycles  simulated,  the 
method  is  not  suggested  for  use  unless  the  expected  cycle  length  is  short 
enough  so  that  an  adequate  number  of  cycles  can  be  simulated  within  one's 
budget  constraint.  If  preliminary  simulation  runs  indicate  that  the 
expected  cycle  length  will  be  excessive  (or  that  the  use  of  the  method 
will  result  in  a variance  addition),  it  is  then  suggested  that  independent 
simulations  be  performed. 
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